import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sympy import *

data = pd.read_csv('data1.csv')
time = np.array(data['time'])
ti = np.array([int(t[:2]) * 60 + int(t[3:]) for t in time])
x = np.array(data['x'])
y = np.array(data['y'])
# 影长
si = np.array([sqrt(xi * xi + yi * yi) for (xi, yi) in zip(x, y)])
# 归一化
si_ = np.array([(s0 + np.mean(si)) / (np.max(si) - np.min(si)) for s0 in si])
st=np.array([(s0 - np.min(si)) / (np.max(si) - np.min(si)) for s0 in si])
# 太阳方位角
omega = np.array([atan(yi / xi) for (xi, yi) in zip(x, y)])
# 太阳直射点纬度
beta = (23.5 * 28 / 91) * 3.1415 / 180
# 北京经度
c2 = 116
# 对h进行经度范围的搜索(步长为1)
h = np.arange(1, 10)
c1 = np.array([np.array(
    [c2 - asin(-sin(o) * cos(t) / cos(beta)) for (o, t) in zip(omega, np.array([atan(hi / x0i) for x0i in s]))]) for hi
    in h])

# 求出c1近似值
c1_ans = np.mean(c1)
Fx = lambda c1, psi, ti: (sin(beta) * sin(psi) - cos(beta) * cos(psi) * cos(
    ti + 3.1415 * 15 * (((c1 - 116) / 15) - 12) / 180)) / (cos(
    asin((sin(beta) * sin(psi) - cos(beta) * cos(psi) * cos(ti + 3.1415 * 15 * (((c1 - 116) / 15) - 12) / 180)))))


omega=np.array([acos(Fx(c1_ans,))])